Última actualización 23 de enero del 2021

Authores: Shalisko V., Castillo-Aja R., Santana E., Valdivia-Ornelas L.

Versión 6.1

Datos fuente

Datos abiertos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127

Fuentes auxiliares

Diccionarios de datos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127 y Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades de INEGI https://www.inegi.org.mx/app/ageeml/

## Warning: package 'zoo' was built under R version 4.0.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sp
## rgdal: version: 1.5-15, (SVN revision 1045)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

ANALISIS MÉXICO

Defunciones por COVID-19 por fecha de inicio de síntomas

ruta_estados <- '../datos_shp/Estados.shp'
estados <- rgdal::readOGR(ruta_estados)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Estados.shp", layer: "Estados"
## with 32 features
## It has 2 fields
ruta_municipios <- '../datos_shp/Municipios.shp'
municipios <- rgdal::readOGR(ruta_municipios)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Municipios.shp", layer: "Municipios"
## with 2456 features
## It has 4 fields
municipios@data$Clave_Mun_Ent_Texto <- as.numeric(as.character(municipios@data$CVE_MUN)) + 
                             1000 * as.numeric(as.character(municipios@data$CVE_ENT))
municipios@data$Clave_Mun_Ent_Texto <- sprintf("%05d", municipios@data$Clave_Mun_Ent_Texto)
#str(municipios@data)
metadatos_m <- datos_municipio[,1:12]
#str(metadatos_m)
dia_m <- datos_municipio[,c(def_fechas)]
dia_antes_m <- cbind(datos_municipio[,c(def_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
incremento_m <- dia_m - dia_antes_m
incremento_m_a100k <- round(100000 * incremento_m / datos_municipio$Pob_Municipio, 2)
incremento_m[is.na(incremento_m)] <- 0 
incremento_m_a100k[is.na(incremento_m_a100k)] <- 0

#dim(incremento_m_a100k)
#head(incremento_m_a100k)
#print(casos_fechas)
#print(casos_fechas_ok)

semanas_lag <- c(seq(from = 0, to = n, by = 7))
def_breaks <- c(0,0.5,1,1.5,2,3,4,5,10,1000000)
def_breaks_labels <- c("0-0.5","0.5-1","1-1.5","1.5-2","2-3","3-4","4-5","5-10",">10")
def_raw_breaks <- c(0,1,2,3,4,5,10,15,20,1000000)
def_raw_breaks_labels <- c("1","2","3","4","5","6-10","11-15","16-20",">20")
def_indice_breaks <- c(-1,-0.75,-0.5,-0.25,-0.1,0.1,0.25,0.5,0.75,1)
def_indice_breaks_labels <- c("menos que esperado (-1)","-0.75","-0.5","-0.25","esperado (0)","0.25","0.5","0.75","mas que esperado (1)")


#for (i in 1:3) {
for (i in 1:46) {
  selection_cols <- paste0("def",fecha_cadena(as.Date(fecha_inicio) + seq(semanas_lag[i]+1, semanas_lag[i]+8)))  
  #print(selection_cols)
  suma_def_semanal <- apply(incremento_m[,selection_cols], 1, sum)
  total_def_semanal <- sum(suma_def_semanal, na.rm = TRUE)
  total_poblacion <- sum(datos_municipio$Pob_Municipio, na.rm = TRUE)
  def_esperados_semanal <- total_def_semanal * datos_municipio$Pob_Municipio / total_poblacion
  ## indice (similar a NDVI)
  def_indice_semanal <- (suma_def_semanal - def_esperados_semanal) / (suma_def_semanal + def_esperados_semanal)

  suma_a100k_semanal <- apply(incremento_m_a100k[,selection_cols], 1, sum)
  datos_semanales <- data.frame(
             Clave_Mun_Ent_Texto = as.character(metadatos_m[,"Clave_Mun_Ent_Texto"]), 
             Municipio = metadatos_m[,"Nom_Mun"],
             Latitud = metadatos_m[,"Lat_Decimal"],
             Longitud = metadatos_m[,"Lon_Decimal"],
             Def_100k = suma_a100k_semanal,
             Def_raw = suma_def_semanal,
             Def_esperado = def_esperados_semanal,
             Def_indice = def_indice_semanal
             )
  
  #print(str(datos_semanales))
  #print(head(datos_semanales, n = 30L))
  mis_municipios <- municipios
  mis_municipios <- merge(municipios, datos_semanales, by = "Clave_Mun_Ent_Texto", all.x = TRUE)
  mis_municipios@data$Class_100k <- as.numeric(cut(mis_municipios@data$Def_100k, 
                                                   breaks = def_breaks))
  mis_municipios@data$Color_100k <- brewer.pal(n = 9, name = 'YlOrRd')[mis_municipios@data$Class_100k]
  mis_municipios@data$Class_raw <- as.numeric(cut(mis_municipios@data$Def_raw, 
                                                   breaks = def_raw_breaks))  
  mis_municipios@data$Color_raw <- brewer.pal(n = 9, name = 'PuRd')[mis_municipios@data$Class_raw]  
  mis_municipios@data$Class_indice <- as.numeric(cut(mis_municipios@data$Def_indice, 
                                                   breaks = def_indice_breaks))  
  mis_municipios@data$Color_indice <- brewer.pal(n = 9, name = 'RdYlBu')[10 - mis_municipios@data$Class_indice]  
  
  #print(head(mis_municipios@data, n = 25L))
  #str(mis_municipios@data)
  
  #str(mis_municipios)
  
  ## grafica defunciones por 100 mil
  plot(mis_municipios,
          col = mis_municipios$Color_100k,
          axes = TRUE, 
          border = NA, 
          #col = "transparent",
          main = paste0("Defunciones semanales por 100 mil habitantes ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright", 
         c("Número de defunciones por COVID-19 con base en los reportes SSA de México"))  
  
  legend("bottomleft", 
         as.character(def_breaks_labels),
         fill = brewer.pal(n = 9, name = 'YlOrRd')[1:9],
         title = "Por 100 mil habitantes"
         )

  ## grafica defunciones
  plot(mis_municipios,
          col = mis_municipios$Color_raw,
          axes = TRUE, 
          border = NA, 
          #col = "transparent",
          main = paste0("Defunciones semanales ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright", 
         c("Número de defunciones por COVID-19 con base en los reportes SSA de México"))  
  legend("bottomleft", 
         as.character(def_raw_breaks_labels),
         fill = brewer.pal(n = 9, name = 'PuRd')[1:9],
         title = "Defunciones semanales"
         )

  ## grafica defuunciones indice
  plot(mis_municipios,
          col = mis_municipios$Color_indice,
          axes = TRUE, 
          border = NA, 
          #col = "transparent",
          main = paste0("Índice de desviación del número de defunciones semanales esperado ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright", 
         c("Áreas en blanco reprezentan zonas sin defunciones registrados en el periodo",
           "Número de defunciones por COVID-19 con base en los reportes SSA de México"))
  legend("bottomleft", 
         as.character(def_indice_breaks_labels),
         fill = brewer.pal(n = 9, name = 'RdYlBu')[9:1],
         title = "Indice de defunciones"
         )  
  
  ## verification plot wit coordinates only
  #plot(x = mis_municipios@data$Longitud,
  #     y = mis_municipios@data$Latitud,
  #     col = mis_municipios@data$Color_raw,
  #     type = "p", pch = 19)
  
}

knitr::asis_output(htmltools::htmlPreserve(cc_html))
Creative Commons License.